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Abstract 

^ I The stability and dynamics of nonlinear Schrodinger superflows past a two- 

dimensional disk are investigated using a specially adapted pseudo-spectral method 
based on mapped Chebychev polynomials. This efficient numerical method allows 
the imposition of both Dirichlet and Neumann boundary conditions at the disk bor- 
!>■ ' der. Small coherence length boundary-layer approximations to stationary solutions 

are obtained analytically. Newton branch-following is used to compute the complete 
bifurcation diagram of stationary solutions. The dependence of the critical Mach 
OO ' number on the coherence length is characterized. Above the critical Mach number, 

. at coherence length larger than fifteen times the diameter of the disk, rarefaction 

O '• pulses are dynamically nucleated, replacing the vortices that are nucleated at small 

coherence length. 
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1 Introduction 



It is well known that, above a critical speed, superfluidity breaks down and dis- 
sipation sets in [1]. Much work has been devoted to the understanding of this 
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phenomenon within the mathematical description of superfluidity provided by 
the nonlinear Schrodinger equation (NLSE) also called the Gross-Pitaevskii 
equation [2,3,4]. 

The NLSE can be used to describe two quite different physical systems: su- 
perfiuid ^He and Bose-Einstein condensates of ultra-cold atomic vapor. 

In the case of superfluid ''He, the NLSE can be considered as a valid math- 
ematical model provided that the temperature is low enough for the normal 
fluid to be negligible. This is clearly the case in recent experiments [5] that 
are performed at temperatures below 130 mK. Note that the excitations of 
superfluid "'He are accurately described by the famous Landau spectrum which 
includes phonons in the low wave number range, and maxons and rotons in 
the high (atomic-scale) wave number range. In contrast, the standard NLSE 
(the equation used in the present article) only contains phonon excitations. 
It therefore incompletely represents the atomic-scale excitations in superfluid 
^He. However, there exist straightforward generalizations of the NLSE [6,7] 
that do reproduce the correct excitation spectrum, at the cost of introducing 
a spatially non-local interaction potential. For reasons of simplicity we shall 
not use such generalizations in the present article. 

Since Bose-Einstein condensation in dilute gases in traps was experimentally 
observed [8,9,10], this field is in rapid evolution: recent results include the pro- 
duction and detection of an isolated quantized vortex [11,12], the nucleation of 
several vortices [13] and details of vortex dynamics [14]. The dynamics of these 
compressible nonlinear quantum fluids is accurately described by the NLSE 
allowing direct quantitative comparison between theory and experiment [15]. 

The stability of Bose-Einstein condensates (BEG) in the presence of a moving 
obstacle can thus be studied in the framework of the NLSE. Raman et al. 
have studied dissipation in a Bose-Einstein condensed gas by moving a blue 
detuned laser beam through the condensate at different velocities [16]. In their 
inhomogcncous condensate, they observed a critical Mach number for the onset 
of dissipation that was compared with the NLSE predictions. 

In their pioneer work, Frisch, Pomeau and Rica [17] performed direct nu- 
merical simulations of the NLSE to study the stability of two-dimensional 
supcrflows around a disk. They observed a transition to a dissipativc regime 
characterized by vortex nucleation that they interpreted in terms of a saddle- 
node bifurcation of the stationary solutions of the NLSE. Later, using nu- 
merical branch-following techniques, Huepe and Brachet [18,19] obtained the 
complete bifurcation diagram in which the stable and unstable branches are 
connected through a saddle-node bifurcation. Asymmetric solutions were also 
found, generated by a secondary pitchfork bifurcation of the stable branch. The 
symmetric and asymmetric unstable solutions correspond respectively to two 
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and one vortices. The critical speed was shown to converge, for small coher- 
ence length, to the Eulerian value computed by Rica [20]. Three-dimensional 
effects leading to a lowering of the critical speed were also considered [21]. 

In all the above numerical studies, the effect of the two-dimensional disk was 
represented in the NLSE by a simple repulsive potential. Thus no boundary 
conditions were applied and the numerical results were (weakly) dependent on 
the details of the repulsive potential. 

One of the main motivations of the present paper is to obtain numerical results 
that are reliable (i.e. do not depend on an ad hoc artificial repulsive poten- 
tial) at finite value of the coherence length. We will thus consider the NLSE 
as a partial differential equation with standard boundary conditions applied 
on the disk. This mathematical problem will be studied by using an efficient 
pseudo-spectral method, based on angular Fourier series and radially mapped 
Chebychev polynomials, that was specifically designed for the present study. 
The numerical solutions will be compared with analytic boundary layer ap- 
proximations, that are valid for small velocity and coherence length. Similar 
expansions were performed for a spherical obstacle in [22] . 

The paper is organized as follows: section 2 contains the governing equations; 
section 3 is devoted to the derivation of the boundary layer analytical ex- 
pressions for Dirichlet conditions; in section 4, we describe the new specially 
designed pseudo-spectral method; section 5 contains validations of the nu- 
merical procedure and new results on bifurcation diagrams and critical Mach 
numbers; in section 6, our results on the dynamically emitted excitations are 
reported, with emphasizing on the nucleation of rarefaction pulses; finally sec- 
tion 7 is our conclusion. More details on the numerical method are found in 
the appendix where the resolutions needed to obtain spectral convergence are 
discussed. 



2 Governing equations 



In this section, we present the hydrodynamic form of the NLSE that models 
the effect of a disk of radius unity (diameter D = 2), moving at constant speed 
V = vBx in a two-dimensional superfluid at rest. In the frame of the disk, the 
system is equivalent to a superflow around a disk, with constant speed — v 
at infinity. Let Q be the plane C deprived of © the disk of radius unity and 
dQ the boundary of the domain, that is the circle of radius unity. We will 
naturally use the polar coordinates (r, 9) such that x = r cos 9 and y = rsm9 
and the associated unit vectors are denoted by {er,eg). The system can then 
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be described with the following action functional 



n 2 



where is a complex field, ^ its conjugate. The speed of sound c and the 
so-called healing length ^ are the physical parameters of the system. is the 
energy of the system that reads 



with 



JQ. 

Ja 2 



- 1) w - (V- - i)V^ 



(2) 

(3) 

(4) 



The presence of the constants —1 in Eq. (4) ensures the convergence of the 
integral. The Euler- Lagrange equation corresponding to (1) provides the NLSE 



idtip = [-e^^ - V' + iv^r^] + iv • , (5) 

defined in the domain Jl. This equation can be mapped into two hydrodynam- 
ical equations by applying Madelung's transformation [1] 

that defines a fluid of density p and velocity 

U = V0 - V (7) 

The real and imaginary parts of the NLSE yield the following equations of 
motion 

9tp + V-(pU) = (8) 

d,<f> = -^(V0)2 + c'il -p) + c'e^ + V • V0. (9) 

These equations correspond respectively to the continuity and the Bernoulli 
equations (with a supplementary quantum pressure term) for a barotropic 

compressible and irrotational flow. Note that two non-dimensional parameters 
control the system: the Mach number Ai = |v|/c (where v is the flow velocity 
at infinity and c the sound speed) and the ratio of the healing length ^ to the 
diameter of the disk D. In the limit ^/D — > 0, the quantum pressure term 
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vanishes and wc recover the system of equations describing an Eulerian flow. 
We now investigate the problem of the boundary conditions on the obstacle. 

In previous studies [19], boundary conditions were applied by adding to the 
NLSE a repulsive potential term strong enough to force the density to zero 
inside the disk. In the present work, we consider the mathematically stan- 
dard Dirichlet and Neumann boundary conditions that will both be directly 
imposed at the border of the obstacle. 

2.1 Dirichlet boundary conditions 

The Dirichlet boundary conditions read V'|r=i = 0. They thus prescribe zero 
density on the obstacle and correspond to the presence of an unpenetrable ob- 
stacle (a laser in a BEC or a solid obstacle in superfluid '^He). They correspond 
to the following conditions, in hydrodynamical variables: first, the condition 
on p is obviously 

p = at r = 1 (10) 

Second, the square root of the density R = y/p being constant on the obstacle, 
we have dtR\r=i = and deR\r=i = 0. The continuity equation (8) expressed 
in term of R then yields drR- U±\r=i = 0, so that the Dirichlet conditions also 
imply 

C/^ = 9^0 - -y cos e = at r = 1 (11) 

2.2 Neumann boundary conditions 

The Neumann boundary, in hydrodynamical variables, read 

drP^O at r = 1 (12) 

U± — drcf) — vcos9 — at r = 1 (13) 

They correspond to the following conditions in term of the complex field ijj: 
driijjexp ( '^^^^'^ ))|r=ro=i — 0. Notc that, compared to the Dirichlet condi- 
tions, the Neumann conditions are more academic than physically realistic. 
Nevertheless, it is interesting to study the influence of such boundary con- 
ditions on the stationary solutions of the problem, especially their effects on 
the boundary layer on the obstacle. For instance, one could think, that with 
such conditions the stationary solution would be "closer" to that of the Eule- 
rian flow than with Dirichlet conditions. We will see below that the situation is 
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more complex. A more physical motivation to study such boundary conditions 
is related to the problem of capillary-gravity free surface flows past a cylin- 
drical obstacle, where the quantum pressure term is replaced by a capillary 
term. In this related problem, the Neumann conditions are physical ones [23]. 



3 Boundary layer solutions — Analytical results 



We now present calculations of the stationary solutions in the limit ^/ D — > 0. 
For non-zero Mach number, 

M = |v|/c (14) 

we define the new phase variable [20] 

— —{(f) — vrcos9)/v. (15) 

The Bernoulli (9) and continuity (8) equations then read 

= e^^-p+l + ^[l-(Vv.)1 (16) 
= pA(/? + Vp • V(^. (17) 

The Dirichlet boundary conditions now read 

p\dQ, = 
dr(fi\dn = 0. 

At finite but small Mach number, we expand p and (p as 

p = + A^2^<i) + . . . + M^fc^W + . . . (18) 
^ = ^(0) + M^ifi^'^ + ■■■ + M^'^ifi^''^ + ■■■ . (19) 

Note that if one knows ip at order Ai^'', on can formally deduce p at order 
^2(fe+i) |-jy solving (16). The potential Lp can then be computed at order 
|-,y solving (17). In order to compute </?, we will have to solve equations 
of the type 

||W + ^^W-^yW = RHS(r) (20) 

Solutions to the corresponding homogeneous equation are 

y{r) ^Ar + Br'^ (21) 
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so that the general equation with nonzero righthandside RHS(r) can be com- 
puted using the method of variation of parameter. Using the boundary con- 
ditions hm^^+oo l/l'") = and dy/dr{r = 1) = yields for the solution of the 
inhomogeneous equation the explicit expression 

y(r) = -— / RHS(ii)(l + u'^)du 

- J /^°^ RHS(M)dM + ^ /^°° v?W^{u)du (22) 

Z Jr 2iT Jr 

provided that the function RHS decreases rapidly enough at infinity. Note that 
the first term of y{r) yields a term of the type C/r. Due to the expressions 
of RHS encountered in the following computations, the two last terms will 
turn out to tend to zero exponentially (on a length scale of order so that 
the behavior at infinity of the function y will be governed by a long-range 
algebraic term that reads 

1 r+oo 

y{r) ~ -tt/ ^^^{u){l + u^)du (23) 

r— >+oo 2r Jl 



We now turn to the computation of the stationary Dirichlet solution. Expres- 
sions for p^^') and are obviously needed to bootstrap the iteration. They 
are obtained by the following considerations. 

When the Mach number is zero, = is solution of the stationary equations 
and p satisfies 

e^-p+1^0 (24) 
Writing p(r, 9) — R'^{r) yields the equation 

C'^AR + R-R^^ ^'^{drr + ^dr)R + R-R^^O (25) 

with boundary conditions -R(l) = 0. A first approximation for the solution of 
this equation, obtained by neglecting the term {^'^/r)drR, reads 

= tanh (^^) (26) 

This result, valid up to order ^, can be improved by setting R = Rq^ + Rf^\ 
Inserting R in (24), collecting the terms of order ^ and solving the resulting 
differential equation yields, after tedious computations, 

^ = [-3 - cosh 2s + (4 + 3s) sech^ s + sinh 2s + 3 tanh s] (27) 
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where s — {r — l)/-\/2^. Thus the exphcit expression 

P^'^=P^^^+P? = {R^^^r + 2RS'^R? (28) 

gives the correct approximation to the density, up to order in the hmit 
^/D ^ 0. 

The velocity potential (^9^°^ satisfies 

A(^(o) = -Vp<°> • V(^<°> + (1 - p<°>)A(^<°>. (29) 

We write (^^°^ = (poller + '^^^^ where v^Euier = (r + 1/r) cos^ is the solution at 
order in M"^ of the Eulerian flow. Using the relation Af^^^jg^, = 0, equation 
(29) yields the following equation for (p^^^ 

A(^<o> = -Vp<°> • vAer - Vp<°> • V(^<°> + (1 - p<°>)A(^<o> (30) 

This equation cannot be solved directly. We thus proceed to a perturbative 

development by writing ip^'^'^ = ipf'^ + (^2°^ where cpf^ is of order ^ and <^2°^ of 
order In the right hand side of equation (30), one can keep at the dominant 
order of our computations the first term and drop the two others. The function 
pf^ is then solution of the equation 

A^? = -VpS°> . wAe. (31) 
The expression of <^f^ can be computed using Eq. (22). Eq. (23) yields [24] 
.(0) 2v/2e-4(log2)e^ 

(fii' ~ — cost' (32) 



r— *+oo 



In order to obtain the full correction at order of the 1/r-algebraic term we 
also need to compute (^2°^ which verifies 

A^f = -Vpf • - • V^f> + (1 - ) A^f (33) 

Using again Eq. (23), a lengthy computation yields 

_^oo^ ^^^^ 



The velocity potential (^^°^ thus reads 



cos^ (35) 



where ip^^l exponentially vanishes at infinity. 
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Note that the compressible Eulerian flow around a disk of radius ri admits at 
order zero in the following solution 

Aer,r,-{r+'-pjcOs9 (36) 

in order to satisfy the boundary condition dr(p\r=ri = 0. Thus, the correction 
to (/^Euiej. is a long-range term that can be physically interpreted as a renor- 
malization of the diameter of the disk : at large distances the superflow is 
equivalent to an Eulerian flow around a disk of radius r^s given by 

The order ^ term was first computed in [24]. Similar results were obtained 
directly, using matched expansions, for a spherical obstacle in [22]. This ref- 
erence also includes the governing matched expansion equations for the case 
of a 2D disk, however the authors did not solve these equations. 

The same procedure with Neumann boundary conditions can be shown to lead 
to a renormalized radius [23] 

(38) 

Note that contrary to the case of Dirichlet conditions, this effective size is de- 
pendent on the Mach number, which was not the case for Dirichlet conditions. 
It is also smaller than the corresponding Dirichlet effective value. 




4 Specially adapted pseudo-spectral method 



We have specifically developed a code that can accurately accommodate both 
large-r asymptotic behavior and thin boundary layers near the obstacle at 
r = 1. It is based on a Chebychev decomposition using an adequate mapping. 
It allows us to consider a unique obstacle in contrast with periodic pseudo- 
spectral methods [19] which in fact model a network of obstacles. 

4-1 Mapping for a unique obstacle 

Using standard polar coordinates {0,r}, together with the relation 

r{z) = z-^ (39) 
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the domain {0 ^ 6 < 2n,—l ^ z ^ 1}, can be mapped into the physical 
domain {x, y}, with + ^ 1. 

The basic mapping is 

X — z'^ cos Q y — z~^ sin 9 (40) 

and the inverse transformation reads 

^ = ^ / '> ^ = arctan(y/a;) + (41) 

Any generic real field ^'(a;, y) (^ stands for Re0, Im'^, etc.) appearing in the 
encountered equations of motion is expressed in the {0,z} domain as 

<if{9,z)^^{x{9,z),y{9,z)) (42) 

with x{9,z) and y{9,z) defined in (40). 

As x{9,z) — x{9 + TT, —z) and y{9,z) — y{9 + tt, —z), the {x,y} domain is 
mapped twice unto the {9,z} domain. A mapped field must therefore satisfy 

^{9,z) = -^{9 + 7r,-z) (43) 



The equations of motion are expressed as partial differential equations in the 
{9, z} domain by writing the differential operators V and A in terms of 9 and 
z derivatives that are polynomial in z, e.g. Aip — z'^^ + z^^^ + ^^ff ■ 



4-2 Generalization of the mapping 



In the special case where ^/D is large, we found useful to generalize the r{z) — 
1/ z mapping to 

rx{zx) = - + (1 - A)^, (44) 

This mapping has the same overall characteristics than 1/z and reduces to 
it for A = 1. For A > 1 it stretches the coordinate, thereby increasing the 
resolution at large distance by moving the collocation points away from the 
r — 1 disk. Generahzations to expressions (40) and (41) are easily derived. 
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4-3 Spatial discretization 



The field ip is spatially discretized, in the {9, z) domain, using a standard 
Chebychev- Fourier pseudo-spectral method [25], based on the expansion 

i^{0,z)^ J2 \j2^n,pT,{z)\ e^pine (45) 

n=l-Ng/2 [p=0 J 

where Tp{z) — cos (p arccos 2;) is the order-p Chebychev polynomial and Nq 
and Nr represent resolutions. 

The pseudo-spectral method calls for using fast Fourier transforms to evaluate 
(45) on the collocation points grid {9m, z^) with 

Q^m<Ne (46) 
Q^k^Nr (47) 

The relation T„(cosx) = cosnx reduces the Chebychev transform appearing in 
(45) to a (fast) Fourier cosine transform. Thus, the evaluation of (45) (and its 
inverse) only requires a time proportional to NgNrlog{NQNr) . Computations 
of nonlinear terms are carried out on the grid representations, while 9 and z 
derivatives are carried out on the Chebychev- Fourier representations. 



9m 


27rm 


Ne ' 




irk 


Zk ■ 





The main virtue of mapping (40) together with expansion (45) is its ability to 
accurately accommodate both large-r asymptotic behavior and thin boundary 
layers near r = 1. Indeed, on the one hand, (45) is an expansion in product 
of polynomials in with functions cosn9 and sinn^^, precisely the type of 
functions needed to capture large-r behavior (see section 5.2 and [20]). On the 
other hand, the accumulation of collocation points Zp (see equation (47)) and 
the regularity of (40) near z — ±1 allows expansion (45) to simultaneously 
resolve boundary layers at r = 1 with thickness of order 1/N^ [25]. 



4-4 Spectral symmetries of the fields 



As ip is real, the coefficients ipn,p in (45) are complex conjugate 

V'-n,p = V'n,p (48) 

They obey an additional relation, stemming from (43). Setting z = cos{9'), 
the fields must be invariant under the transformation 9^^9 + 11, 9' 9' + 7r. 
In spectral space, this transformation reads ipn,p ^ (~l)"(~l)^^n,p; implying 
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(49) 



Thus the ■ipn,p coefficients are non-zero only when {n,p) are jointly even or 
jointly odd. This relation, similar to that found in the Taylor- Green vortex 
[26], is used to speed-up the evaluation of (45) by a factor 2, using specially 
designed even-odd Fast Fourier transforms. 

Integral of mapped fields are performed on the collocation points using the 
discrete formula 

j^rdrd9^{r,9) = - — —^ ^ i^{9m, z,).Ji^^-{zM^p) (50) 



5 Stationary solutions — Numerical results 

This section is devoted to the numerical determination of stationary solutions 
using the branch-following method detailed in the appendix. We first focus on 
the particular case of the Eulerian flow (that is when ^/D — 0). This case has 
been previously investigated using methods based on series in Mach number 
by Rica [20], and the critical Mach number is known with great precision. 
We next compare analytical results of section 3 with numerically obtained 
profiles of boundary layers with Dirichlet conditions. It is thus a good test of 
the numerical precision and efficiency of our new method, presented above in 
section 4. 

The rest of the section contains the numerical results on the bifurcation dia- 
grams and the stationary solutions of the NLSE at small and large coherence 
lengths, for the two types of boundary conditions: Dirichlet and Neumann. 
We discuss the dependence on ^/D of the critical Mach number. 

5.1 Eulerian limit 

In the limit ^/D — > 0, the NLSE turns into the equations of an Eulerian 
compressible flow 

9*0 = ~{Vcf>f + c'{l - p) + V • V0 (51) 
dtp = -pA(f) - Vp • V(/) + V • Vp (52) 

that are respectively the Bernoulli and continuity equations. We now search 
for their stationary solutions. Note that, knowing the stationary fleld 0, the 
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Bernoulli equation yields explicitly an expression of p that reads 

P=l-^(V0)' + ^v-V0 (53) 

Therefore, is solution to the following equation 

= -pA0 - Vp • V0 + V • Vp (54) 

with density p given by Eq. (53) and a unique boundary condition on the disk 
(instead of two in the NLSE case) 

U±_ — dr(f) — V COS ^ = at r =1 (55) 

Using the branch following method presented in the appendix yields the nu- 
merical stationary solutions of the two-dimensional Eulerian flow with respect 
to the Mach number. The critical Mach number is then the one at which the 
local Mach number 



|U| _ |V0 
ci°^ y/p 



A^'~=^ = ^^- (56) 



reaches 1, at {x — 0,y — ±1)[20]. 

The value of the computed critical Mach number determined in this way de- 
pends on the resolution. It is found to decrease when the resolutions in 9 
and r increase. In order to obtain 11 significant digits, the (minimum) needed 
resolution is Nq x Nr = 512 x 32. The critical Mach number then found is 
Aic = 0.36969705259(9). In order to obtain the same precision as that of se- 
ries methods [20] {A4f^'^^ = 0.36969(7)), it is sufficient to use the resolution 
X Nr — 128 X 16, that is only 8 radial grid points in physical space. 



5.2 Comparison with analytical boundary layer results for Dirichlet condtions 



We now compare the analytical results of section 3 with numerically obtained 
profiles of boundary layers with Dirichlet conditions. Figure 1(a) displays 
boundary layer profiles of the density square-root {R — ^/p) computed at 
C/D = 1/200 and M ^ 0. 

To stress the agreement between analytical and numerical results, it is more 
convenient to substract the term (26) i?o°^ — tanh((r — l)/^/2^) in the nu- 
merical profiles and compare the higher order terms thus obtained with the 
analytical expression. 

In the same way. Figure 1(b) presents the order variation of the effective 
radius as the following combination 6ef£ = ((^eff/ro)^ — l)/(^/ro)). The line 
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2\/2+ [10 — 16 log 2) /3] (^/ ro) is shown on the same graph since expression (37) 
predicts that 



3eff 



^ ^ 10-161og2/^\ ^,^2N 
2\/2+ — ^1 ^ I +0(e^) 



(57) 



The value of {res/roY was extracted from the numerical results by calculat- 
ing the coefficient in cos9/r of the velocity potential of the stationary state 
substracted by the corresponding Eulerian flow coefficient (see Eq. (36)). 

The agreement between analytical and numerical results is very good for small 
^/D emphasizing the ability of our method to compute thin boundary layers. 




(b) 



0.02 0.04 0.06 0.08 



Fig. 1. (a) Plot of Rrisir)) with s{r) = (r - 1)/V2^ (see Eq. (27) and below) 
together with the function R™^^{r) — tanh((r — l)/-\/2^) where i?°'^™(r) is the nu- 
merical result obtained with our numerical method for ^/D = 1/200 and M = 0. 
The agreement is excellent, (b) Calculation of 6 function of ^ together with 

the curve 2\/2 + (10 — 161og 2)(^/ro)/3 (see text). The difference between the two 
curves is due to the term in SeS of higher order in ^. Note that the agreement is 
very good for small ^/rg. 



5.3 Bifurcation diagrams and stationary states at small coherence length 



We present the bifurcation diagrams and the stationary solutions of the NLSE 
at small coherence length, for the two types of boundary conditions: Dirichlet 
and Neumann. 

The numerical methods presented in section 4 and the appendix also converge 
very well in the NLSE case. For instance for Dirichlet conditions, the resolution 
needed to compute a whole bifurcation diagram is lower than in previous 
studies by Huepe and Brachet [19]. With the present method, the resolution 
needed in the case = 1/20 is A^g x A^^ = 64 x 64 whereas Huepe et al. 
needed a spatial (rectangular) resolution x Ny = 256 x 128 for the same 
ratio ^/D. The gain in resolution is then of a factor 8. This factor increases 
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for smaller ^/D. These excellent convergence properties are detailed in the 
appendix, section A. 4. 2. 



5.3.1 Bifurcation diagrams 

In order to study bifurcation diagrams, we define a new free energy by: 

^[^,^]=^o[^,^]-v-y2cc/ d£n^(^-^), (58) 

with n = — the unit vector normal to the boundary. The rightmost term 
in (58) does not affect the equation of motion and is always zero for the 
■01 an — Dirichlet boundary. For Neumann conditions, this term ensures that 
a stationary solution i/jq is an actual extremum of the functional JF, i.e., T 
satisfies J?^['?/'o + '^'0? "^o + ^"0] ~ -^["^Oi "^o] = at first order in b'll). This property 
implies the existence of a generic cusp in T at the bifurcation point (see 
figure 2). 

For simplicity, we will use the notation T{M) = J-'[ipo{M.),ipo{Ai)]. The 
values of JF(A1) — JF(0) (the change of energy JF, relative to zero Mach number) 
is displayed in figure 2 as a function of the Mach number Ai for various values 
oi^/D and the two types of boundary conditions. As can be seen by inspection 
of the figure, for each ^/-D, the stable branch (solid fine) disappears with the 
unstable solution (dashed line) at a saddle- node bifurcation when Ai = Ale- 
There are no stationary solutions beyond this point. This qualitative behavior 
is the signature of a Hamiltonian saddle node bifurcation. 

By inspection of figure 2(b), we can see that the stable stationary branches 
for Neumann conditions are almost superimposed on the Eulcr branch which 
is not the case for Dirichlet conditions. This is due to the fact that Dirichlet 
conditions impose a zero of the density at the border of the disk, contrary to 
Neumann conditions and Eulerian flow. 

In Fig. 2(a), at a Mach number smaller than A4c, the unstable symmetric 
branch (dashed line, circle, ^/D = 1/20) bifurcates at a pitchfork to a pair of 
asymmetric branches (dotted line, C,/D = 1/20) [18]. It can be directly checked 
on our results (see figure 3, middle) that the secondary pitchfork bifurcation 
breaks the y i— > —y symmetry of the flow for both boundary conditions. 

At a fixed Mach number, the energy difference between a stable and an un- 
stable solution corresponds to the energy barrier necessary to dynamically 
nucleate an excitation. Note that this barrier for a symmetric unstable solu- 
tion is about twice that of an asymmetric unstable solution. 
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Fig. 2. Bifurcation diagrams for small coherence lengths. Energy functional 
J^{M) — J-'{0) versus Mach number, (a) Dirichlet conditions; (b) Neumann con- 
ditions. For S^/D = 1/20, the asymmetric unstable solution branch is represented (it 
stands for a one- vortex branch stemming from a pitchfork bifurcation). At a fixed 
Mach number, the energy difference between a stable and an asymmetric unstable 
solution is roughly half the energy difference between a stable and a symmetric 
unstable solution. 

5.3.2 Stationary solutions 

By visualizing the stationary solutions of the NLSE, the branches of Fig. 2 can 
be related to the presence of vortices. Figure 3 shows the density p = IV'P of 
typical stationary solutions for M. = 0.3 and C^/D = 1/20 for the two types of 
boundary conditions. It is apparent by inspection of the figure that the stable 
branch is irrotational (figure 3, top) while the asymmetric unstable branch 
corresponds to a one-vortex solution (figure 3, middle) and the symmetric 
unstable branch, to a two- vortex solution (figure 3, bottom). 

For Dirichlet boundary conditions, similar results were found with periodic 
pseudo-spectral codes [19]. However, our method directly imposes the cor- 
rect boundary conditions without resorting to an artificial repulsive potential. 
Also note that the critical Mach number is here determined for a single obsta- 
cle, whereas a periodic array of obstacles was used in previous study. Huepe 
et al. [19] find for ratio ^/D = 1/40, A1^"<=p<= ~ 0.3817 whereas we obtain 
Aic — 0.3941. A single obstacle perturbs less the fiow than an infinite ar- 
ray of obstacles (even with large separation), it is therefore natural to find a 
higher critical Mach number in our simulations. As M. is increased, the dis- 
tance between the vortices and the obstacle for the unstable branches (figures 
3, middle, bottom, Dirichlet) decreases. At a certain Mnv < A^c, the vortices 
disappear on the surface on the cylinder, generating an irrotational fiow (see 
[18] for a detailed study of the Mach number at which one or two vortices 
emerge from the disk) . 

Note that the branch following procedure used to compute the unstable branches 
bifurcating from the stable branch conserves the velocity circulation. The total 
velocity circulation around the disk is null. The two- vortex solution conserves 
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Dirichlet stable 



Neumann stable 




Dirichlet 1 vortex Neumann 1 vortex 




Dirichlet 2 vortices Neumann 2 vortices 




Fig. 3. Density p = of stationary solutions for $^/D = 1/20 and M = 0.3 far 
from the bifurcation threshold: (top) stable solution, (middle) asymmetric unstable 
solution and (bottom) symmetric unstable solution. Left : Dirichlet conditions; right: 
Neumann conditions. 



the total circulation since the two vortices are counter-rotating. For the one- 
vortex solution, an image vortex located at the middle of the obstacle has to 
be invoked. This point will be reconsidered in section 5.4. 
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5.3.3 Variation of critical Mach number with ^/D 

We now study the dependence on ^/D of the critical Mach number J\4c- 

Our numerical method needs a slight modification to allow us to explore the 
large $,/D regime. The transformation r{z) is modified such that mesh-points 
situated near the obstacle are stretched (sec section 4, equation (44)). This 
procedure avoids wasting resolution close to the cylinder. 
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Fig. 4. Critical Mach number M.^ versus ^jD. Note that Dirichlet boundary con- 
dition solutions admit a smaller A^c than Neumann boundary condition solutions. 
They both tend to the Euler critical Mach number as ^/D decreases. 

Results are displayed in figure 4. For a given type of boundary conditions, 
the Mach number decreases with ^/D and, for both boundary conditions, it 
converges towards the Euler hmit for small ^/D. 

First note that the value of the critical Mach number is lower than 1. As we 
are interested in stationary solutions with density approaching 1 at infinity 
like polynomials in 1/r (see sections 3 and 4) , the speed of the obstacle v has 
to remain below the speed of sound c. Otherwise radiation of sound waves 
would occur in the same way as discussed in [27,28] . 

We now discuss the case of Dirichlet boundary conditions and will extend the 
argument to Neumann conditions at the end of this section. 

For Dirichlet conditions, we have seen in section 3 and in section 5.2 that, 
in the case of small ^/D, the boundary layer has a thickness of order ^. The 
situation at large ^/i? is quite different. 

We now show that the effect of the cylinder on the fiow at r > 1 is vanishingly 
small when ^ — > oo. At zero Mach number, the density p^{r) = -R|(r) satisfies 
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^2 (^drrR^ + ^drR^^ - i?| + i?^ = (59) 

with -R^(l) = and R^{+oo) — 1. For large r/^, = — 1 satisfies after 
linearization of (59) 

(drrR^ + ^drR^^ - 2R^ = (60) 

Asymptotically for large r/^, we have R^{r) ~ 1 + ^|pp''°^ with 



Ko(4^) 

papprox _ ^ ^ 



for a given constant /j,^. 

Using a shooting method, we have numerically solved Eq. (59) starting from 
r = [B is a sufficiently large constant) with initial conditions i?'^"™(i?^) = 
1 + ^^PP'"°"(50 and 9^i?™'^(5^) - a^^|PP™"(S^) and adjusting the constant 
/^^ so that i?""''^(l) = 0. 

Our shooting method indicates that 

lim /i£ = 1+ (62) 
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Fig. 5. Plot of the function i?^ ~ 1 + ji^pp^°^ together with i?^"™ both expressed in 
the variable s = (r — 1)/^ for S^/D = 2500 (in this case, ~ 1.057). The function 
1 + ]i^PP^°^ is a very good approximation of the solution of equation (59) except 
close to the cylinder (s = 0). 

Figure 5 displays the function 1 + ^Ipp'""'' together with the numerical solu- 
tion of Eq. (59) calculated by the shooting method, expressed in term of the 
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rescaled variable s — (r — l)/^. Thus, at large ^/D and for s > 0, 1 + is 
a very good approximation of the solution of Eq. (59). Points close to the cylin- 
der differ, which is obvious since 1 + ^|pp'°^ does not vanish at the cylinder. 
Therefore we have for fixed s > 

pn^^^ 1 Ko(f + V2s) Ko{V2s) 

i?t(l + ts) - 1 — Aif ^ — 1= ~ 7^ — T-^ (63) 

^ ^ Ko(f ) logC ^ ' 

In this sense, the square-root of the density approaches logarithmically 
the uniform state. Then the effects of the obstacle on the flow at large ^/D 
for Dirichlet conditions are very small. Thus the critical velocity is expected 
to increase with ^/D and the larger ^/D, the closer to 1 the critical Mach 
number. 

As the Neumann conditions perturb even less the fluid than the Dirichlet con- 
ditions, one can easily understand why the critical Mach number for Neumann 
conditions increases faster at large ^/D than for Dirichlet conditions. 

Turning now to the small ^/D regime, the critical Mach number for Neumann 
conditions is also found to be larger than that of Dirichlet conditions. This 
point is quite surprizing since, at small Mach number, stationary solutions 
for Neumann conditions approach the Euler stationary states better than the 
stationary solutions for Dirichlet conditions (see the bifurcation diagrams on 
figure 2). However, we can offer the following semi-quantitative argument. 
The critical Mach number decreases with decreasing $,/D. We have shown in 
section 3 that the effective radius rcfj(^) of stationary NLS flows at small Mach 
numbers was bigger for Dirichlet conditions than for Neumann conditions. 
Assuming that this result holds for bigger Mach numbers (of order A4^^^^^), 
one can consider that the Neumann conditions stationary solutions have the 
same critical Mach number as the Dirichlet stationary solutions when they 
reach the same ratio ^/Des{^), imposing therefore smaller values oi ^/D in 
the Neumann case. 

Finally, note that we have found no numerical indication showing that 
A^c(Dirichlet) could become bigger that A4c(Neumann), for very small values 
of ^/D. 



5.4 Bifurcation diagrams and stationary states at large coherence length 



The large ^/D regime could be reached experimentally by considering BEC 
with large coherence lengths perturbed by a sharply focused detuned laser. 
As seen in the previous section, the critical Mach number tends very quickly 
towards 1 for Neumann boundary conditions. Furthermore, these conditions 
are academic and have no experimental equivalent in BEC. Thus, we will 
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study the limit of large ^/D only for Dirichlet conditions, the experimentally 
realistic ones. 

The bifurcation diagram, computed for Dirichlet boundary conditions and 
different ^/D, is displayed on figure 6(a). Just like in the small ^/D case, 
a branch of stable solutions is connected to a branch of unstable solutions 
through a saddle-node bifurcation at a critical Mach number M.^- The values of 
A^c are seen to approach 1 when ^/ D increases. The corresponding stationary 
solutions are also displayed on figure 6. The stable solution (b) is irrotational 
while the two unstable solutions contain respectively one (d) or two vortices (c) 
far from the critical Mach number. However, close enough to the bifurcation 
tip, the unstable stationary solution shows no 27r phase jump and therefore 
no vortices are present. We will come back to this point in the next section. 

As already pointed in section 5.3.2, the one-vortcx solution that breaks the 
symmetry y ^ —y associates a vortex located outside the obstacle to an 
image vortex situated inside the obstacle. The image vortex is clearly visible 
on figure 6(d) because its core is larger than the obstacle itself (compare 
with figure 6b or c). Note that the large ^/D energy difference between the 
stable and the asymmetric solution branches is not half the energy difference 
between the stable and the symmetric solution branches, contrary to the small 
^/D case. This stems from the fact that, for large ^/D, the energy of the 
asymmetric branch also includes the additional contribution of the now visible 
image vortex. 



6 Dynamical results 

Solutions of the NLSE {in the absence of an obstacle) in dimension 2, moving 
at constant speed while preserving their shape, have been exhibited by Roberts 
et al. [29,30]. These solutions are pairs of counterrotating vortices but also 
what they called rarefaction pulse (depletion pulse with non zero density and 
therefore no vorticity). 

A natural question is then to know which kind of excitations can be nucleated 
past a disk. 

6. 1 Nucleation of vortices 

The stationary solutions obtained numerically provide us with adequate initial 
data for the study of dynamical solutions. Indeed, after a small perturbation, 
their integration in time will generate a dynamical evolution with very small 
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Fig. 6. (a) Bifurcation diagram for large ^jD for Dirichlet conditions. Density of 
stationary solutions for ^jD = 20 and ^A = 0.25 far from the bifurcation threshold 
{M.C — 0.86): (b) stable solution, (c) symmetric unstable solution and (d) asym- 
metric unstable solution. 

acoustic emission. This procedure also provides an efficient way to start vor- 
tical dynamics in a controlled manner. 

It is already known from studies performed using a repulsive potential (see 
figure 6 of [19]) that, at small values of ^/D, vortex pairs are dynamically 
nucleated. The same behavior is obtained using the present numerical method 
with Dirichlet boundary conditions (data not shown). This behavior persists 
when using Neumann boundary conditions, as shown on figure 7 that displays 
the nucleation of a + and — vortex pair (which will be followed by a periodic 
emission of other pairs). The phase of the complex field exhibits a 27r jump 
for each vortex (see figure 7(d)). 

6.2 Nucleation of rarefaction pulse 

For large ^/D, using Dirichlet boundary conditions, we have proceeded in the 
same way as in the previous subsection by perturbing an unstable symmetric 
solution at Mach number M > Aic to observe the nature of the nucleated 
excitations. The behavior is somewhat more complicated. We show on Table 1 
the nature of the emitted excitations as a function oi ^/D and J^/M.^- Foi' 
^/D > 15, and obstacle speed above Ale, a rarefaction pulse (RP) is dynam- 
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Fig. 7. Typical vortex nuclcation in the case of Neumann boundary conditions, for 
^/D = 1/20. (a,b,c) Plot of the density at different times (in arbitrary unit t = 0, 
t = 230, t = 260 respectively), (d) At t = 260, phase of the system. Note the 
27r-phase jumps around the two points where the density vanishes: iso-phase lines 
emerge from these two points. 

ically obtained; for ^/D < 15, there exists a threshold in M. under which 
vortex pair rather than a rarefaction pulse emerges. In some cases (BL for 
border line), it is unclear whether we have an emitted rarefaction pulse or 
a vortex pair (the density minimum in such limit cases approaches zero and 
there is a strong variation of the phase). By measuring the borderline speed 
of translation v (in the frame at rest) of the emitted excitation (see Table 2), 
we found that it is very close to the known limit speed of translation -uy at 
which occurs the change in the nature of excitations in 2D superflow in the 
absence of an obstacle: v^/c = 0A3\/2 ~ 0.61 [30]. 

A desexcitation of an unstable stationary solution at $,/D = 17.5 creates a 
rarefaction pulse as shown in figure 8. The minimum of density of such a pulse 
is non zero (here the minimum equals approximately 0.081), and no phase 
jump is present (see figure 8(d)). 

The change in the nature of the excitation can be understood by the following 
qualitative argument. A rarefaction pulse can be seen as the superposition of 
a pair of vortices so close to each other that the minimum density is non zero. 
Close to the critical Mach number, no vortex is detached from the disk. When 
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nucleated, these vortices follow the boundary of the obstacle and then leave it 
separated by a distance of order D. For large ^, the vortices are so close that 
they give rise to a rarefaction pulse. 

We did not study the periodic emission of rarefaction pulses at supercritical 

regime as done in previous studies [19]. Our numerical method is not adapted 
for such problems: the mesh is more and more stretched far from the obstacle so 
that one lacks resolution at long distance to resolve the nucleated excitations. 
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Tabic 1 

Phase diagram of the nature of emitted excitations as a function of the ratio ^/D and 
the Mach number normalized by the critical Mach number. VP and RP respectively 
stand for vortex pair and rarefaction pulse. BL stands for limit cases where it is 
hard to distinguish the exact nature of the excitation (the density minimum is very 
close to zero and the phase has a strong variation). 
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Tabic 2 

Speed of translation v/c oi the nucleated excitations at M/ Mc = 1.1 as a function 
of ^/D. For the limit case ^/D = 10, we find v ~ 0.59; the change in the nature 
of excitation in a 2D superflow without obstacle appears at speed equal to fy/c ~ 
0.61 [30]. 



7 Conclusion 

The main virtue of the pseudo-spectral method that we have used in the 
present study is its ability to accurately accommodate both large-r asymp- 
totic behavior and thin boundary layers near the cylindrical obstacle, at r = 1. 
Indeed, using modest resolutions, we were able to obtain the Eulerian critical 
Mach number with 11 significant digits. In the NLSE case, spectral conver- 
gence was obtained on the whole bifurcation diagram for values oi^/ D as low 
as 1/120. 
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(a) 



(b) 




Fig. 8. Nucleation of a rarefaction pulse with Dirichlet boundary conditions for 
^/D = 17.5 and M./M.C = 1-1: density at different times (in arbitrary unit) (a) 
f = 0, (b) t = 150, (c) t = 220. (d) Plot of the phase at t = 220. There is no phase 
jump, hence no vorticity. The speed of translation of this rarefaction pulse in frame 
at rest is 0.80 and its minimum of density is approximately 0.081. 

Small coherence length boundary-layer approximations to the stationary solu- 
tions were calculated. These analytical results were found to be in very good 
agreement with the numerical results. The long-range contribution was phys- 
ically interpreted as a renormalization of the diameter of the disk. 

As a by product of our new method, we were able to investigate not only the 
physically realistic case of Dirichlet boundary conditions but also the more 
academic case of Neuman conditions. The influence of the boundary condi- 
tions on the stationary solutions of the problem, especially their effects on the 
boundary layer and on the critical Mach number, were investigated. 

For Dirichlet boundary conditions the qualitative results previously obtained, 
using periodic pseudo-spectral codes [19], were recovered. However, our new 
method directly imposes the correct boundary conditions, without resorting to 
an artificial repulsive potential. Also, the newly obtained critical Mach number 
is here determined for a single obstacle, whereas a periodic array of obstacles 
was considered in previous studies. Thus, the present article presents the first 
precise quantitative determination of the critical Mach number as a function 
of ^/D in this reference problem. 
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Finally wc were able to show that a transition occurs in the nature of the 
emitted excitation at large coherence length. For ^/D > 15, the nucleated 
excitations are rarefaction pulses whereas, aX ^/D < 15, both vortices and 
rarefaction pulses can be obtained. 
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A Numerical Methods 

We detail here the numerical procedures used in the simulations together 
with the numerical convergences with spatial resolutions for the critical Mach 
number and the stationary solutions. 

A.l Spectra 

In this article, we define the r-spectrum and ^-spectrum of a field ip spectrally 
represented by ipn,p as the respective sequence of numbers 

Ng 

SP'(P)= E |^n,pr O^p^Nr (A.l) 

Sp^(n) = \^n,p\' O^n^^ (A.2) 

Only half the ^-spectrum is considered for the ^-representation is complex- 
conjugated. 



A.2 Implementation oj the boundary conditions 

Contrary to the analytical computations, we work with the complex variable 
ip in order to take account of the possible presence of vortices. In order to 
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impose the boundary conditions, we set 

Tn cos 9 



V2c^ 



zr 



(A.3) 
(A.4) 



Note that this change of variables does not affect the Dirichlet conditions that 
still read 



With these new variables, Neumann conditions read 
and the NLSE turns into 



(A.5) 



(A.6) 



V2^ 



-e^AV'm + (|V'mr-l)^/'m + iv • Vz/^, 



with ip- 



m |r=l 



+ ^ [e'(V$o)'V'm + e'2i(V$o) VV^^] - V • (V$o)V'm (A.7) 
(Dirichlet) or drip,^\r=i = (Neumann). 



A.3 Time steppings 



A.3.1 Stationary States 

We search for stationary solutions of the dynamics equations (5) or Eq. (51- 
52). Note that stationary solutions are those of the equivalent diffusive equa- 
tions that read in the abbreviated form 



— = L* + W(*) 
at 

In the general case (^ 7^ 0) , we have 



V2^ 



L = A 



(A.8) 

(A.9) 
(A.IO) 



f (V$o)'V^n + f 2i(V$o)VV^„ + V • (V$o)V'n 



In the Eulerian case ({ = 0), as discussed in section 5.1, these definitions 
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reduce to 



* = L = A W = V(^V0) - V • V(/) (A.ll) 

with 

^=-^(V</>)' + v-V0 (A.12) 

To integrate (A. 8) a mixed implicit-explicit first-order time stepping scheme 
is used: 

*(i + r) = (I-rL)-i(I + TW)*(i) (A.13) 
where I is the identity operator and r the time step. 

The Helmholtz operator (1— rL), block-diagonal with respect to Fourier modes, 
is easily inverted in the Fourier-Chebychev representation using the LU algo- 
rithm [31]. 

As called for the r method [25] the boundary conditions (55), (A. 5) or (A. 6) 
are substituted to the equations (A. 8) for the highest Chebychev modes Ti^^_i 
and Tjv^. The operator (I — rL) is thus modified before inversion. 

This relaxation method can only reach stable stationary solutions of (A. 8). 
In order to also capture unstable stationary solutions [32] we use the Newton 
branch- following method detailed in [33], [18], [19]. 



A. 3. 2 Branch following procedure 

We search for fixed points of (A.13), a condition strictly equivalent to the 
stationarity of (A. 8). Each Newton step requires solving a hnear system for 
the decrement to be subtracted from ^: 



(I-tL)-^ (I + tDW)-I]V 

(I-TL)-^(I + rW) - I 



(A. 14) 



where DW(^) is the Frechet derivative, or Jacobian matrix, of W evaluated 
at . Equation (A. 14) is equivalent to: 



(I - tL)-V(L + DW)V' = (I - tL)-V(L + W)* 



(A.15) 



The role of r is formally that of the time step in (A.13), but in (A. 14) or 
(A.15), its value can be taken to be arbitrarily large. For r — > oo, (A.15) 
becomes: 



L-^(L + DW)V' = L-^(L + W)* 



(A.16) 
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In order to solve the linearized systems stemming from the Newton method, we 
use BiCGSTAB [34]. We vary r empirically to optimize the preconditioning 
and convergence of BiCGSTAB. A few hundred BiCGSTAB iterations are 
usually required to solve the linear system. 

A. 3. 3 Dynamics 

We write Eq. (5) in the abbreviated form: 

^ = L'* + W'(*) (A.17) 
at 

where 

^ = ^^ L' = -iL w' = -iW(^^) (A.18) 

with W defined in equation (A. 10) 

Equation (A.17) is time stepped using the implicit Euler scheme 

= (1 - tL)-' + tW„] (A.19) 

The boundary conditions (55), (A. 5) or (A. 6) are imposed by modifying the op- 
erator (I— tL)~^, as done for the relaxation time stepping algorithm (A. 8) [25]. 

A. 4 Numerical convergence 
A. 4.1 Euler 

Table A.l shows the error on our reference Mach number versus and A^^- 
Note that the errors are mainly due to a lack of Fourier modes in 9. Thus, 
when a sufficient number of Fourier modes is reached, increasing the resolution 
in r yields a better precision. 
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Table A.l 



Relative error versus resolution on the critical Mach number calculated by taking 
as a reference Mc = 0.36969705259(9) calculated at (512 x 32). 
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Due to the use of Chebychev polynomials, the boundary layers of the NLS flow 
computed at low Mach number are well resolved thanks to the large number of 
collocation points at the vicinity of the obstacle, as we will see in next section. 



A. 4-2 Small coherence length solutions 

Our numerical method based on Chebychev polynomial expansions allows to 
solve the boundary layer of order ^ by refining the collocation points near the 
boundary conditions. The smaller ^, the larger the radial resolution Nr must 
be. The azimuthal resolution No depends also on the value of ^ through the 
multiplication of complex fields with a phase term such that ^o{9, r) = v 

and ip-ca = '^e^*°. The phase term $o is inversely proportional to ^ and needs 
sufficient No points in order to be resolved. Table A. 2 lists the resolutions used 
for computing the bifurcation diagram for each ^/D. Spectral convergence is 
achieved for all stationary solutions as shown in figure A.l. 
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Table A.2 

Azimuthal and radial resolutions used for computing the bifurcation diagram for 
different ^/D for the two types of boundary conditions. 
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Fig. A.l. Stationary solution spectra with = 1/20, A^i x A^ = 64 x 64, and 
Neumann boundary conditions: (a) spectra and (b) r spectra for a stable solution 
far from the bifurcation, for the solution at the bifurcation, and for a one vortex and 
a two vortex unstable solutions. Spectral convergence is achieved for all stationary 
solutions. 



A. 4- 3 Large coherence length solutions 

For large ratio $,/D, we have modified the mapping in order to stretch in the 
radial direction the collocation points next to the obstacle (see equation (44)). 
The choice of the value of A depends on the value oi ^/ D and the resolution 
Nq X Nr of the system. 
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Figure A. 2 shows azimuthal and radial spectra for a symmetric unstable sta- 
tionary solution at ^/D = 20 with a dilatation parameter A = 80 for two 
different resolutions. Spectral convergence is achieved for all stationary solu- 
tions. 




Fig. A. 2. Symmetric unstable stationary solution spectra with ^/D = 20, for two 

resolutions {Ng x Nr = 128 x 128 et 128 x 512) and Dirichlet boundary conditions: 
(a) 6 spectra and (b) r spectra. Spectral convergence is achieved for all stationary 
solutions. 
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